## Replication script for "The Political Economy of Bureaucratic Overload: Evidence from Rural Development Officials in India"
## Aditya Dasgupta (adasgupta3@ucmerced.edu) and Devesh Kapur (dkapur1@jhu.edu)

## Notes: All required files attached in enclosed folder. See readme.txt file for description of data sets and variable names.
## Once appropriate working directory is assigned, this R script should run in its entirety to produce all of the main figures
## and tables in the paper. All required packages listed below and may need to be installed. Please run script in its entirety
## as some of the downstream/later analyses depends on objects produced during earlier analyses. 

## Working directory
####################

rm(list=ls())
setwd("C:/Users/Adi/Dropbox/Devesh_BDO_project/Revise_Resubmit/final_documents/replication_final_files") ## Set to working directory with enclosed data

## Libraries
############

library(ggplot2)
library(ggridges)
library(DiagrammeR)
library(gam)
library(lfe)
library(msm)
library(xtable)
library(stargazer)
library(rsvg)
library(DiagrammeRsvg)
library(tm)
library(stargazer)
library(lfe)
library(rdrobust)
library(ivpack)
library(wordcloud)

IKbandwidth <-function (X,Y,cutpoint=NULL,verbose=FALSE,kernel="triangular") {
  #Implementation of Imbens-Kalyanaraman optimal bandwidth
  # for regression discontinuity
  sub<-complete.cases(X)&complete.cases(Y)
  X <- X[sub]
  Y <- Y[sub]
  Nx<-length(X)
  Ny<-length(Y)
  if(Nx!=Ny)
    stop("Running and outcome variable must be of equal length")
  if(is.null(cutpoint)) {
    cutpoint<-0
    if(verbose) cat("Using default cutpoint of zero.\n")
  } else {
    if(! (typeof(cutpoint) %in% c("integer","double")))
      stop("Cutpoint must be of a numeric type")
  }
  #Now we should be ready to start
  #Pilot bandwidth
  h1<-1.84*sd(X)*Nx^(-1/5)
  left<-X>=(cutpoint-h1) & X<=cutpoint
  right<-X>cutpoint & X<=(cutpoint+h1)
  Nl<-sum(left)
  Nr<-sum(right)
  Ybarl<-mean(Y[left])
  Ybarr<-mean(Y[right])
  fbarx<-(Nl+Nr)/(2*Nx*h1)
  varY<-(sum((Y[left]-Ybarl)^2)+sum((Y[right]-Ybarr)^2))/(Nl+Nr)
  medXl<-median(X[X<=cutpoint])
  medXr<-median(X[X>cutpoint])
  Nl<-sum(X<cutpoint)
  Nr<-sum(X>=cutpoint)
  cX<-X-cutpoint
  if(sum(X[left]>medXl)==0 | sum(X[right]<medXr)==0)
    stop("Insufficient data in vicinity of the cutpoint to calculate bandwidth.")
  #Model a cubic within the pilot bandwidth
  mod<-lm(Y~I(X>=cutpoint)+poly(cX,3,raw=T),subset=(X>=medXl&X<=medXr))
  m3<-6*coef(mod)[5]
  #New bandwidth estimate
  h2l<-3.56*(Nl^(-1/7))*(varY/(fbarx*max(m3^2,0.01)))^(1/7)
  h2r<-3.56*(Nr^(-1/7))*(varY/(fbarx*max(m3^2,0.01)))^(1/7)
  left<-(X>=(cutpoint-h2l)) & (X<cutpoint)
  right<-(X>=cutpoint) & (X<= (cutpoint+h2r))
  Nl<-sum(left)
  Nr<-sum(right)
  if(Nl==0 | Nr==0)
    stop("Insufficient data in vicinity of the cutpoint to calculate bandwidth.")
  #Estimate quadratics for curvature estimation
  mod<-lm(Y~poly(cX,2,raw=T),subset=right)
  m2r<-2*coef(mod)[3]
  mod<-lm(Y~poly(cX,2,raw=T),subset=left)
  m2l<-2*coef(mod)[3]
  rl<-720*varY/(Nl*(h2l^4))
  rr<-720*varY/(Nr*(h2r^4))
  #Which kernel are we using?
  # Method for finding these available in I--K p. 6
  if(kernel=="triangular") {
    ck<-3.43754
  } else if (kernel=="rectangular") {
    ck<-5.40384
  } else if(kernel=="epanechnikov") {
    ck<-3.1999
  } else if(kernel=="quartic" | kernel=="biweight") {
    ck<-3.65362
  } else if(kernel=="triweight") {
    ck<-4.06065
  } else if(kernel=="tricube") {
    ck<-3.68765
  } else if(kernel=="gaussian") {
    ck<-1.25864
  } else if(kernel=="cosine") {
    ck<-3.25869
  } else {
    stop("Unrecognized kernel.") 
  }
  #And there's our optimal bandwidth
  optbw<-ck*(2*varY/(fbarx*((m2r-m2l)^2+rr+rl)))^(1/5)*(Nx^(-1/5))
  left<-(X>=(cutpoint-optbw)) & (X<cutpoint)
  right<-(X>=cutpoint) & (X<= (cutpoint+optbw))
  if(sum(left)==0 | sum(right)==0)
    stop("Insufficient data in the calculated bandwidth.")
  names(optbw)<-NULL
  if(verbose) cat("Imbens-Kalyanamaran Optimal Bandwidth: ",sprintf("%.3f",optbw),"\n")
  return(optbw)
}

########################################################
##                  						##
##                  Figures and Tables			##
##									##
########################################################

## Figure 1: Theoretical Diagram
################################

DiagrammeR::grViz("digraph {

graph [layout = dot]

# define the global styles of the nodes. We can override these in box if we wish
node [shape = rectangle, style = filled, fillcolor = Linen]

process1 [label =  'Unclear \n Responsibility', fillcolor = Beige]
process2 [label =  'Weak Electoral Incentives to \n Invest in Local Bureaucracy', fillcolor = Beige]
process3 [label = 'Under-resourced and \n Overloaded Bureaucracy']
process4 [label = 'Inability to \n Specialize']
results [label= 'Poor Implementation']

# edge definitions with the node IDs
process1 ->process2 -> process3 ->results
process3 -> process4 -> results
}")  %>%
    export_svg %>% charToRaw %>% rsvg_pdf("fig1.pdf") ## Output diagram

## Figure 2: Rural Development Blocks in Survey
###############################################

## Please contact authors for survey coordinates as this involves identifiable information

## Figure 3: Focus group data: Discussion on Resources
######################################################

focus_group <- read.csv("focus_group.csv") ## Load focus group data

pdf(file="fig3.pdf", width=10, height=7)

par(mfrow=c(1,2))

responses <- 100*table(focus_group$resources_score)/sum(table(focus_group$resources_score)) ## Barplot of quantitative scores
names(responses) <- c("Score 1", "Score 2", "Score 3", "Score 4", "Score 5")
barplot(responses, xlab="Focus Group Resources Score", ylab="Percentage")

docs <- c(as.character(focus_group$resources_notes)) ## create vector of documents
docs <- Corpus(VectorSource(docs)) ## create corpus from vector of documents
docs <- tm_map(docs, stripWhitespace) ## no white space
docs <- tm_map(docs, tolower) ## lower case
docs <- tm_map(docs,removePunctuation) ## remove punctuation
docs <- tm_map(docs, removeWords, stopwords("english")) ## remove stop words
docs <- tm_map(docs, stemDocument) ## stem words
dtm <- DocumentTermMatrix(docs) ## create term-doc matrix
dtm_sparse <- removeSparseTerms(dtm, .95)## up to 95% missing allowed
TDM1 <- as.matrix(dtm_sparse) #Convert this into a matrix format
frequency1 <- sort(colSums(TDM1), decreasing = TRUE) #Gives you the frequencies for every word
set.seed(94611)## set seed
wordcloud(names(frequency1), frequency1, col="blue",random.order=FALSE,scale=c(4,1))

dev.off()

## Table 1: Bureaucratic Resource Index 
#######################################

tab1 <- read.csv("tab1.csv")

names <- c("Full-time Employees", "Contract Employees",
	"Vehicles", "Computers")
means <- c(mean(tab1$fullstaff), mean(tab1$contractstaff),
	mean(tab1$vehicles), mean(tab1$computers))
sds <- c(sd(tab1$fullstaff), sd(tab1$contractstaff),
	sd(tab1$vehicles), sd(tab1$computers))
lower <- c(quantile(tab1$fullstaff,p=.25), quantile(tab1$contractstaff, p=.25),
	quantile(tab1$vehicles, p=.25), quantile(tab1$computers, p=.25))
upper <- c(quantile(tab1$fullstaff,p=.75), quantile(tab1$contractstaff, p=.75),
	quantile(tab1$vehicles, p=.75), quantile(tab1$computers, p=.75))
loadings <- c(fit$loadings[1,1],fit$loadings[2,1],fit$loadings[3,1],
	fit$loadings[4,1])
mat <- cbind(means,sds,lower, upper, loadings)
write.table(print(xtable(mat), 
	"html", include.rownames=FALSE, 
               html.table.attributes='align="left"'), file="table1.htm", row.names=F, col.names=F, quote=F)

## Figure 4: Patterns of Daily Time Allocation by Rural Development Officials
#############################################################################

holder <- read.csv("fig4.csv")

## Create ridgeline plot

jpeg(file="fig4.jpg",width=800,height=800)

par(mar=c(5,7,4,1)+.1)
plot(holder[,7]+1.4, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n", xaxt="n", ylab="", xlab="Time",
	xlim=c(1.7,20.3))
polygon(c(1:21,21:1), c(holder[,7]+1.4, c(rep(1.4,times=21))),col="grey")
axis(1, labels=c("9 am",10,11,12,1,2,3,4,5,6,"7 pm"),
	at=c(1,3,5,7,9,11,13,15,17,19,21))
axis(2, labels=c("Other","Citizens","Planning","Field","Management","Politicians","Forms","Unrelated Duties"),
	at=c(0,.2,.4,.6,.8,1,1.2,1.4), las=2)
abline(h=c(0,.2,.4,.6,.8,1,1.2,1.4), lty=3, col="lightgrey")
abline(v=c(1,3,5,7,9,11,13,15,17,19,21), lty=3, col="lightgrey")

points(holder[,1]+1.2, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,1]+1.2, c(rep(1.2,times=21))),col="grey")

points(holder[,6]+1, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,6]+1, c(rep(1,times=21))),col="grey")

points(holder[,2]+.8, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,2]+.8, c(rep(.8,times=21))),col="grey")

points(holder[,4]+.6, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,4]+.6, c(rep(.6,times=21))),col="grey")

points(holder[,3]+.4, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,3]+.4, c(rep(.4,times=21))),col="grey")

points(holder[,5]+.2, type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,5]+.2, c(rep(.2,times=21))),col="grey")

points(holder[,8], type='l', lwd=3, ylim=c(0,1.6), col=1, yaxt="n")
polygon(c(1:21,21:1), c(holder[,8], c(rep(0,times=21))),col="grey")

dev.off()

## Figure 5: Differences in Block-level Clarity of Responsibility
#################################################################

## Please contact authors for survey coordinates as this involves identifiable information

## Table 2: Impact on under-resourcing on the performance of programs
######################################################################

tab2 <- read.csv("tab2.csv")

## Days of Employment

## Cross-state

lm1 <- felm(nrega ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | India | 0 | State.Name,
		data=tab2)

## Within-state

lm2 <- felm(nrega ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		data=tab2)

## Within-district

lm3 <- felm(nrega ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		data=tab2)

## Expenditures

## Cross-state

lm4 <- felm(expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | India | 0 | State.Name,
		data=tab2)

## Within-state

lm5 <- felm(expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		data=tab2)

## Within-district

lm6 <- felm(expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		data=tab2)

## Expenditures

## Cross-state

lm7 <- felm(labor_expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | India | 0 | State.Name,
		data=tab2)

## Within-state

lm8 <- felm(labor_expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		data=tab2)

## Within-district

lm9 <- felm(labor_expenditures ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		data=tab2)

write.table(stargazer(lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm8,lm9, omit=c("District.Name","State.Name"), type="html"), 
	file="table2.htm", row.names=F, col.names=F, quote=F) ## Output Table 2


## Table 3: Impact on a non managerially intensive program (Swachh Bharat)
##########################################################################

tab3 <- read.csv("tab3.csv")

## Number of toilets per 1000 capita

## Cross-state

lm1 <- felm(toilets ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | India | 0 | State.Name,
		data=tab3)

## Within-state

lm2 <- felm(toilets ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		data=tab3)

## Within-district

lm3 <- felm(toilets ~ normresources + 
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		data=tab3)

write.table(stargazer(lm1,lm2,lm3, omit=c("District.Name","State.Name"), type="html"), 
	file="table3.htm", row.names=F, col.names=F, quote=F) ## Output Table 3

## Table 4: Impact of Resources on Bureaucratic Behavior
########################################################

tab4 <- read.csv("tab4.csv")

## Cross-state

lm1 <- felm(I(forms*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight, data=tab4)

lm2 <- felm(I(management*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm3 <- felm(I(planning*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm4 <- felm(I(field*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm5 <- felm(I(politicians*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm6 <- felm(I(citizens*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm7 <- felm(I(unrelated*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm8 <- felm(I(forms*100 + management*100 + planning*100) ~ normresources + corruption +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm9 <- felm(ELF ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | 0 | 0 | State.Name,
		weights=tab4$weight,data=tab4)

write.table(stargazer(lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm9, omit=c("District.Name","State.Name"), type="html"), 
	file="table4a.htm", row.names=F, col.names=F, quote=F) ## Output Table 4a

## Within-state

lm1 <- felm(I(forms*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm2 <- felm(I(management*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm3 <- felm(I(planning*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm4 <- felm(I(field*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm5 <- felm(I(politicians*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm6 <- felm(I(citizens*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm7 <- felm(I(unrelated*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm8 <- felm(I(forms*100 + management*100 + planning*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm9 <- felm(ELF ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | State.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

write.table(stargazer(lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm9, omit=c("District.Name","State.Name"), type="html"), 
	file="table4b.htm", row.names=F, col.names=F, quote=F) ## Output Table 4b

## Within-district

lm1 <- felm(I(forms*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm2 <- felm(I(management*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm3 <- felm(I(planning*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm4 <- felm(I(field*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm5 <- felm(I(politicians*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm6 <- felm(I(citizens*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm7 <- felm(I(unrelated*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm8 <- felm(I(forms*100 + management*100 + planning*100) ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

lm9 <- felm(ELF ~ normresources +
	joinpromotion + temporary + TOTAL + MINORITY + male + BA + postgrad + literacy + hav_distance | District.Name | 0 | State.Name,
		weights=tab4$weight,data=tab4)

write.table(stargazer(lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm9, omit=c("District.Name","State.Name"), type="html"), 
	file="table4c.htm", row.names=F, col.names=F, quote=F) ## Output Table 4c

## Table 5: Determinants of Bureaucratic Resources
##################################################

tab5 <- read.csv("tab5.csv")

## Panel A: Across States

lm1 <- felm(normresources ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance | 0 | 0 | State.Name,
		data=tab5)

lm2 <- felm(fullstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance | 0 | 0 | State.Name,
		data=tab5)

lm3 <- felm(contractstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | 0 | 0 | State.Name,
		data=tab5)

lm4 <- felm(vehicles ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | 0 | 0 | State.Name,
		data=tab5)

lm5 <- felm(computers ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | 0 | 0 | State.Name,
		data=tab5)

write.table(stargazer(lm1, lm2, lm3, lm4, lm5, omit=c("District.Name","State.Name"), type="html"), 
	file="table5a.htm", row.names=F, col.names=F, quote=F) ## Output Table 5a

## Panel B: Within states

lm1 <- felm(normresources ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | State.Name | 0 | State.Name,
		data=tab5)

lm2 <- felm(fullstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance | State.Name | 0 | State.Name,
		data=tab5)

lm3 <- felm(contractstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | State.Name | 0 | State.Name,
		data=tab5)

lm4 <- felm(vehicles ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | State.Name | 0 | State.Name,
		data=tab5)

lm5 <- felm(computers ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | State.Name | 0 | State.Name,
		data=tab5)

write.table(stargazer(lm1, lm2, lm3, lm4, lm5, omit=c("District.Name","State.Name"), type="html"), 
	file="table5b.htm", row.names=F, col.names=F, quote=F) ## Output Table 5b

## Panel C: Within districts

lm1 <- felm(normresources ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance   | District.Name | 0 | State.Name,
		data=tab5)

lm2 <- felm(fullstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance | District.Name | 0 | State.Name,
		data=tab5)

lm3 <- felm(contractstaff ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance   | District.Name | 0 | State.Name,
		data=tab5)

lm4 <- felm(vehicles ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | District.Name | 0 | State.Name,
		data=tab5)

lm5 <- felm(computers ~ GIS_share_pop_aligned + I(1-GIS_legislator_population_fractionalization) + MINORITY + hav_distance  | District.Name | 0 | State.Name,
		data=tab5)

write.table(stargazer(lm1, lm2, lm3, lm4, lm5, omit=c("District.Name","State.Name"), type="html"), 
	file="table5c.htm", row.names=F, col.names=F, quote=F) ## Output Table 5c

## Figure 6: Close Elections Regression Discontinuity Design
############################################################

fig6 <- read.csv("fig6.csv")

jpeg(file="fig6a.jpg",width=420,height=420)

rdrobust::rdplot(fig6$GIS_share_pop_aligned, fig6$ruling_margin, c = 0, p = 3, 
	y.lim=c(0,1), nbins=5,
		x.lim=c(min(fig6$ruling_margin,na.rm=T),max(fig6$ruling_margin,na.rm=T)), 
			title="", x.label="Ruling Party Win/Loss Vote Share Margin in Constituency Containing Block HQ",
				y.label="Aligned Legislator (Share Population)", lines.lwd=3)
dev.off()

jpeg(file="fig6b.jpg",width=420,height=420)

rdplot(fig6$normresources-mean(fig6$normresources), fig6$ruling_margin, c = 0, p = 3,
	y.lim=c(-.5,.5),
		x.lim=c(min(fig6$ruling_margin,na.rm=T),max(fig6$ruling_margin,na.rm=T)), 
			title="", x.label="Ruling Party Win/Loss Vote Share Margin in Constituency Containing Block HQ", 
				y.label="Bureaucratic Resource Index", cex=5, lines.lwd=3, nbins=5)
dev.off()

## Table 6: Fuzzy RDD Estimates
###############################

tab6 <- read.csv("tab6.csv")
tab6$State.Name <- as.character(tab6$State.Name) ## Make sure loaded as character vector and not factor

## Panel A: Cubic Polynomial

iv1 <- ivreg(normresources ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab6)

iv2 <- ivreg(fullstaff ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab6)

iv3 <- ivreg(contractstaff ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab6)

iv4 <- ivreg(vehicles ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab6)

iv5 <- ivreg(computers ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab6)

write.table(stargazer(cluster.robust.se(iv1, tab6$State.Name), cluster.robust.se(iv2, tab6$State.Name),
	cluster.robust.se(iv3, tab6$State.Name), cluster.robust.se(iv4, tab6$State.Name), 
		cluster.robust.se(iv5, tab6$State.Name), digits=2, type="html"), 
	file="table6a.htm", row.names=F, col.names=F, quote=F) ## Output Table 6a

## Panel B: Local-linear regression

IK <- IKbandwidth(tab6$ruling_margin, tab6$normresources, kernel="triangular")

iv1 <- ivreg(normresources ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right , data = tab6[abs(tab6$ruling_margin)<= IK,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= IK])

iv2 <- ivreg(fullstaff ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab6[abs(tab6$ruling_margin)<= IK,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= IK])

iv3 <- ivreg(contractstaff ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab6[abs(tab6$ruling_margin)<= IK,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= IK])

iv4 <- ivreg(vehicles ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab6[abs(tab6$ruling_margin)<= IK,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= IK])

iv5 <- ivreg(computers~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab6[abs(tab6$ruling_margin)<= IK,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= IK])

write.table(stargazer(cluster.robust.se(iv1, tab6$State.Name[abs(tab6$ruling_margin)<= IK]), cluster.robust.se(iv2, tab6$State.Name[abs(tab6$ruling_margin)<= IK]),
	cluster.robust.se(iv3, tab6$State.Name[abs(tab6$ruling_margin)<= IK]), 
		cluster.robust.se(iv4, tab6$State.Name[abs(tab6$ruling_margin)<= IK]), 
			cluster.robust.se(iv5, tab6$State.Name[abs(tab6$ruling_margin)<= IK]), digits=2, type="html"), 
	file="table6b.htm", row.names=F, col.names=F, quote=F) ## Output Table 6b

## Panel C: Local-polynomial robust estimates with triangular kernel

iv1  <- rdrobust(tab6$normresources, tab6$ruling_margin, cluster=tab6$State.Name, 
	fuzzy=tab6$GIS_share_pop_aligned, p=1, q=2, kernel="tri")
bandwidth1 <- iv1$bws[1,]
bandwidth2 <- iv1$bws[2,]

iv2 <- rdrobust(tab6$fullstaff, tab6$ruling_margin, cluster=tab6$State.Name, 
		fuzzy=tab6$GIS_share_pop_aligned, bwselect="cerrd", p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv3 <- rdrobust(tab6$contractstaff, tab6$ruling_margin, cluster=tab6$State.Name, 
		fuzzy=tab6$GIS_share_pop_aligned, bwselect="cerrd", p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv4 <- rdrobust(tab6$vehicles, tab6$ruling_margin, cluster=tab6$State.Name, 
		fuzzy=tab6$GIS_share_pop_aligned, bwselect="cerrd", p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv5 <- rdrobust(tab6$computers, tab6$ruling_margin, cluster=tab6$State.Name, 
		fuzzy=tab6$GIS_share_pop_aligned, bwselect="cerrd", p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

Coefs <- t(matrix(c(round(iv1$coef, digits=2)[3],round(iv2$coef, digits=2)[3],round(iv3$coef, digits=2)[3],
	round(iv4$coef, digits=2)[3],round(iv5$coef, digits=2)[3])))
SEs <- t(matrix(c(round(iv1$se, digits=2)[3],round(iv2$se, digits=2)[3],round(iv3$se, digits=2)[3],
	round(iv4$se, digits=2)[3],round(iv5$se, digits=2)[3])))

write.table(stargazer(rbind(Coefs, paste("(",SEs,")",sep="")), type="html"), 
	file="table6c.htm", row.names=F, col.names=F, quote=F) ## Output Table 6c

## Panel D: Local difference in means

iv1 <- ivreg(normresources ~ GIS_share_pop_aligned | 
	aligned, data = tab6[abs(tab6$ruling_margin)<= bandwidth1[1]/2,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= bandwidth1[1]/2])

iv2 <- ivreg(fullstaff ~ GIS_share_pop_aligned | 
	aligned, data = tab6[abs(tab6$ruling_margin)<= bandwidth1[1]/2,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= bandwidth1[1]/2])

iv3 <- ivreg(contractstaff ~ GIS_share_pop_aligned | 
	aligned, data = tab6[abs(tab6$ruling_margin)<= bandwidth1[1]/2,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= bandwidth1[1]/2])

iv4 <- ivreg(vehicles ~ GIS_share_pop_aligned | 
	aligned, data = tab6[abs(tab6$ruling_margin)<= bandwidth1[1]/2,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= bandwidth1[1]/2])

iv5 <- ivreg(computers ~ GIS_share_pop_aligned | 
	aligned, data = tab6[abs(tab6$ruling_margin)<= bandwidth1[1]/2,], weights=tab6$IK_weight[abs(tab6$ruling_margin)<= bandwidth1[1]/2])

write.table(stargazer(cluster.robust.se(iv1, tab6$State.Name[abs(tab6$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv2, tab6$State.Name[abs(tab6$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv3, tab6$State.Name[abs(tab6$ruling_margin)<= bandwidth1[1]/2]),
			cluster.robust.se(iv4, tab6$State.Name[abs(tab6$ruling_margin)<= bandwidth1[1]/2]),
				cluster.robust.se(iv5, tab6$State.Name[abs(tab6$ruling_margin)<= bandwidth1[1]/2]), digits=2, type="html"), 
	file="table6d.htm", row.names=F, col.names=F, quote=F) ## Output Table 6d

## Table 7: Fuzzy RDD Estimates of Balance on Covariates
########################################################

tab7 <- read.csv("tab7.csv")
tab7b <- read.csv("tab7b.csv") ## One less observation due to missing literacy data
tab7$State.Name <- as.character(tab7$State.Name) ## Make sure loaded as character and not factor vector
tab7b$State.Name <- as.character(tab7b$State.Name) ## Make sure loaded as character and not factor vector

## Panel A: Cubic polynomial

iv1 <- ivreg(I(TOTAL/1000) ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv2 <- ivreg(MINORITY ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv3 <- ivreg(hav_distance ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv4 <- ivreg(literacy ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7b)

iv5 <- ivreg(N ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv6 <- ivreg(I(1-GIS_legislator_population_fractionalization) ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv7 <- ivreg(male ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv8 <- ivreg(joinpromotion ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv9 <- ivreg(postgrad ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv10 <- ivreg(temporary ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv11 <- ivreg(age ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv12 <- ivreg(years ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv13 <- ivreg(promotion ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv14 <- ivreg(empowered ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

iv15 <- ivreg(power ~ GIS_share_pop_aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right| 
	aligned + polynomial1 + polynomial2 + polynomial3 + polynomial1_right + polynomial2_right + polynomial3_right,  data = tab7)

write.table(stargazer(cluster.robust.se(iv1, tab7$State.Name), cluster.robust.se(iv2, tab7$State.Name),
	cluster.robust.se(iv3, tab7$State.Name), cluster.robust.se(iv4, tab7b$State.Name), 
		cluster.robust.se(iv5, tab7$State.Name),
	cluster.robust.se(iv6, tab7$State.Name), cluster.robust.se(iv7, tab7$State.Name),
		cluster.robust.se(iv8, tab7$State.Name), cluster.robust.se(iv9, tab7$State.Name), 
			cluster.robust.se(iv10, tab7$State.Name), 
				cluster.robust.se(iv11, tab7$State.Name),
	cluster.robust.se(iv12, tab7$State.Name), cluster.robust.se(iv13, tab7$State.Name),
		cluster.robust.se(iv14, tab7$State.Name), cluster.robust.se(iv15, tab7$State.Name), digits=2, type="html"), 
	file="table7a.htm", row.names=F, col.names=F, quote=F) ## Output Table 7a

## Panel B: Local linear regression approach

iv1 <- ivreg(I(TOTAL/1000) ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv2 <- ivreg(MINORITY ~  GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv3 <- ivreg(hav_distance ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv4 <- ivreg(literacy ~  GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7b[abs(tab7b$ruling_margin)<= IK,], weights=tab7b$IK_weight[abs(tab7b$ruling_margin)<= IK])

iv5 <- ivreg(N ~  GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv6 <- ivreg(I(1-GIS_legislator_population_fractionalization) ~  GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv7 <- ivreg(male ~  GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv8 <- ivreg(joinpromotion ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv9 <- ivreg(postgrad ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv10 <- ivreg(temporary ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv11 <- ivreg(age ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv12 <- ivreg(years ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv13 <- ivreg(promotion ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv14 <- ivreg(empowered ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

iv15 <- ivreg(power ~ GIS_share_pop_aligned + running_left + running_right | 
	aligned + running_left + running_right, data = tab7[abs(tab7$ruling_margin)<= IK,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= IK])

write.table(stargazer(cluster.robust.se(iv1, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), cluster.robust.se(iv2, tab7$State.Name[abs(tab7$ruling_margin)<= IK]),
	cluster.robust.se(iv3, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), cluster.robust.se(iv4, tab7b$State.Name[abs(tab7b$ruling_margin)<= IK]), 
		cluster.robust.se(iv5, tab7$State.Name[abs(tab7$ruling_margin)<= IK]),
	cluster.robust.se(iv6, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), cluster.robust.se(iv7, tab7$State.Name[abs(tab7$ruling_margin)<= IK]),
		cluster.robust.se(iv8, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), cluster.robust.se(iv9, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), 
			cluster.robust.se(iv10, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), 
				cluster.robust.se(iv11, tab7$State.Name[abs(tab7$ruling_margin)<= IK]),
	cluster.robust.se(iv12, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), cluster.robust.se(iv13, tab7$State.Name[abs(tab7$ruling_margin)<= IK]),
		cluster.robust.se(iv14, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), 
			cluster.robust.se(iv15, tab7$State.Name[abs(tab7$ruling_margin)<= IK]), digits=2, type="html"), 
	file="table7b.htm", row.names=F, col.names=F, quote=F) ## Output Table 7b

## Panel C: Robust local linear estimates

iv1  <- rdrobust(tab7$TOTAL/1000, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv2  <- rdrobust(tab7$MINORITY, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv3  <- rdrobust(tab7$hav_distance, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv4  <- rdrobust(tab7b$literacy, tab7b$ruling_margin, cluster=tab7b$State.Name, 
	fuzzy=tab7b$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv5  <- rdrobust(tab7$N, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv6  <- rdrobust(1-tab7$GIS_legislator_population_fractionalization, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv7  <- rdrobust(tab7$male, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv8  <- rdrobust(tab7$joinpromotion, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv9  <- rdrobust(tab7$postgrad, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv10  <- rdrobust(tab7$temporary, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv11  <- rdrobust(tab7$age, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv12  <- rdrobust(tab7$years, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv13  <- rdrobust(tab7$promotion, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv14  <- rdrobust(tab7$empowered, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

iv15  <- rdrobust(tab7$power, tab7$ruling_margin, cluster=tab7$State.Name, 
	fuzzy=tab7$GIS_share_pop_aligned, p=1, q=2, kernel="tri",
			h=bandwidth1, b=bandwidth2)

Coefs <- t(matrix(c(round(iv1$coef, digits=2)[3],round(iv2$coef, digits=2)[3],round(iv3$coef, digits=2)[3],
	round(iv4$coef, digits=2)[3],round(iv5$coef, digits=2)[3],
		round(iv6$coef, digits=2)[3], round(iv7$coef, digits=2)[3], round(iv8$coef, digits=2)[3], 
			round(iv9$coef, digits=2)[3], round(iv10$coef, digits=2)[3],round(iv11$coef, digits=2)[3],
	round(iv12$coef, digits=2)[3], round(iv13$coef, digits=2)[3], 
			round(iv14$coef, digits=2)[3], round(iv15$coef, digits=2)[3])))
SEs <- t(matrix(c(round(iv1$se, digits=2)[3],round(iv2$se, digits=2)[3],round(iv3$se, digits=2)[3],
	round(iv4$se, digits=2)[3],round(iv5$se, digits=2)[3],
		round(iv6$se, digits=2)[3], round(iv7$se, digits=2)[3], round(iv8$se, digits=2)[3], 
			round(iv9$se, digits=2)[3], round(iv10$se, digits=2)[3], round(iv11$se, digits=2)[3],
	round(iv12$se, digits=2)[3], round(iv13$se, digits=2)[3], round(iv14$se, digits=2)[3],
		round(iv15$se, digits=2)[3])))

write.table(stargazer(rbind(Coefs, paste("(",SEs,")",sep="")), type="html"), 
	file="table7c.htm", row.names=F, col.names=F, quote=F) ## Output Table 7c

## Panel D: Local Difference in means

iv1 <- ivreg(I(TOTAL/1000) ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv2 <- ivreg(MINORITY ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv3 <- ivreg(hav_distance ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv4 <- ivreg(literacy ~ GIS_share_pop_aligned | 
	aligned, data = tab7b[abs(tab7b$ruling_margin)<= bandwidth1[1]/2,], weights=tab7b$IK_weight[abs(tab7b$ruling_margin)<= bandwidth1[1]/2])

iv5 <- ivreg(N ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv6 <- ivreg(I(1-GIS_legislator_population_fractionalization) ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv7 <- ivreg(male ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv8 <- ivreg(joinpromotion  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv9 <- ivreg(postgrad  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv10 <- ivreg(temporary  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv11 <- ivreg(age  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv12 <- ivreg(years  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv13 <- ivreg(promotion  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv14 <- ivreg(empowered  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

iv15 <- ivreg(power  ~ GIS_share_pop_aligned | 
	aligned, data = tab7[abs(tab7$ruling_margin)<= bandwidth1[1]/2,], weights=tab7$IK_weight[abs(tab7$ruling_margin)<= bandwidth1[1]/2])

write.table(stargazer(cluster.robust.se(iv1, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv2, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv3, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
cluster.robust.se(iv4, tab7b$State.Name[abs(tab7b$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv5, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv6, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
cluster.robust.se(iv7, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv8, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv9, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
cluster.robust.se(iv10, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv11, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv12, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
cluster.robust.se(iv13, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
	cluster.robust.se(iv14, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]),
		cluster.robust.se(iv15, tab7$State.Name[abs(tab7$ruling_margin)<= bandwidth1[1]/2]), digits=2, type="html"), 
	file="table7d.htm", row.names=F, col.names=F, quote=F) ## Output Table 7d


########################################################
##                  						##
##                  End of Script				##
##									##
########################################################






